rm(list = ls())
library(tidyverse)
library(broom)
library(haven)
library(polycor)
library(wCorr)
library(quanteda)
library(quanteda.textplots)
require(glmnet)
require(quanteda.textstats)
library(ggcorrplot)
library(tidymodels)

load("../working/het_effects_stability.Rdata")
load("../working/het_effects_constraint.Rdata")
load("../working/het_effects_polarization.Rdata")

# Plot issue-specific effects for attention

pol_plot <- polarization_het_effects$attention %>%
  bind_rows(.id = "level") %>%
  mutate(name = factor(name, levels = name[level == "High"][order(est[level=="High"])])) %>%
  mutate(level = case_when(level == "High" ~ "High Attention",
                           level == "Low" ~ "Low Attention"),
         level = factor(level, levels = c("Low Attention", "High Attention"))) %>%
  ggplot(aes(x = est, xmin = low, xmax = high, y = name, col = level)) + 
  geom_pointrange(position = position_dodge(.5)) + 
  theme_bw() + 
  geom_vline(xintercept = 0, linetype = 2) + 
  scale_color_manual("", values = c("gray", "black")) + 
  xlab("Effect of Reason Giving") + 
  ylab("") + ggtitle("Polarization")+ 
  theme(legend.position = "bottom")

stab_plot <- stability_het_effects$attention %>%
  bind_rows(.id = "level") %>%
  mutate(name = case_when(name == "MinimumWage" ~ "Minimum Wage",
                          name == "Offensive" ~ "Offensive Speech",
                          name == "Tax" ~ "High Tax",
                          name == "Transgender" ~ "Trans Rights",
                          name == "Unemployment" ~ "Unemployment Support",
                          name == "ZeroHours" ~ "Zero Hours")) %>%
  mutate(name = factor(name, levels = rev(c("Unemployment Support", "Minimum Wage", "High Tax", "Trans Rights", "Zero Hours", "Offensive Speech")))) %>%
  mutate(level = case_when(level == "High" ~ "High Attention",
                           level == "Low" ~ "Low Attention"),
         level = factor(level, levels = c("Low Attention", "High Attention"))) %>%
  ggplot(aes(x = est, xmin = low, xmax = high, y = name, col = level)) + 
  geom_pointrange(position = position_dodge(.5)) + 
  theme_bw() + 
  geom_vline(xintercept = 0, linetype = 2) + 
  scale_color_manual("", values = c("gray", "black")) + 
  xlab("Effect of Reason Giving") + 
  ylab("") + ggtitle("Stability") + 
  theme(legend.position = "bottom")


const_plot <- constraint_het_effects$attention %>%
  bind_rows(.id = "level") %>%
  mutate(name = factor(name, levels = name[level == "High"][order(est[level=="High"])])) %>% 
  mutate(level = case_when(level == "High" ~ "High Attention",
                           level == "Low" ~ "Low Attention"),
         level = factor(level, levels = c("Low Attention", "High Attention"))) %>%
  ggplot(aes(x = est, xmin = low, xmax = high, y = name, col = level)) + 
  geom_pointrange(position = position_dodge(.5)) + 
  theme_bw() + 
  geom_vline(xintercept = 0, linetype = 2) + 
  scale_color_manual("", values = c("gray", "black")) + 
  xlab("Effect of Reason Giving") + 
  ylab("") + ggtitle("Constraint") + 
  theme(legend.position = "bottom")

pdf("../out/outcomes/heterogeneous_effects_attention.pdf", width = 15, height = 5)
gridExtra::grid.arrange(const_plot, pol_plot, stab_plot, ncol = 3)
dev.off()

png("../out/outcomes/heterogeneous_effects_attention.png", width = 9000, height = 3000, res = 600)
gridExtra::grid.arrange(const_plot, pol_plot, stab_plot, ncol = 3)
dev.off()

# Mean effects for all covariates

get_mean_effects <- function(effect_list){
  lapply(1:length(effect_list), function(i) 
    lapply(effect_list[[i]], function(x) data.frame(mean = unique(x$Diff_mean),
                                                    low = unique(x$Diff_low),
                                                    high = unique(x$Diff_high))) %>%
      bind_rows() %>%
      mutate(var = names(effect_list)[[i]],
             names = names(effect_list[[i]]))) %>%
    bind_rows()
}

polarization_het_effects <- get_mean_effects(polarization_het_effects)
stability_het_effects <- get_mean_effects(stability_het_effects)
constraint_het_effects <- get_mean_effects(constraint_het_effects)

polarization_het_effects$outcome <- "Polarization"
stability_het_effects$outcome <- "Stability"
constraint_het_effects$outcome <- "Constraint"

het_effects <- bind_rows(polarization_het_effects, stability_het_effects, constraint_het_effects)
het_effects$y <- 1:nrow(het_effects)

het_effects <- tibble(het_effects) %>%
  mutate(var = case_when(var == "attention" ~ "Attention",
                         var == "gender" ~ "Gender",
                         var == "age" ~ "Age",
                         var == "education" ~ "Education",
                         var == "eu_vote" ~ "Brexit",
                         var == "vote19" ~ "Vote19"),
         names = ifelse(names == "Other" & var == "Brexit", "Did not vote", names),
         var = factor(var, levels = c("Attention", "Age", "Gender", "Education", "Brexit","Vote19")),
         names = tidytext::reorder_within(names, y, var))



ggplot(het_effects, aes(x = mean, y = names, xmin = low, xmax = high)) + 
  geom_pointrange() + 
  facet_grid(var ~ outcome, scales = "free_y") + 
  theme(strip.background = element_blank(), strip.placement = "outside") + 
  theme_bw() + 
  geom_vline(xintercept = 0, lty = 2) + 
  xlab("Conditional Average Treatment Effect") + 
  ylab("") + tidytext::scale_y_reordered()

ggsave("../out/outcomes/heterogeneous_effects.pdf", width = 11, height = 5.5)


# Plot average effects for treatment duration heteogeneity

load(file = "../working/dur_effects_constraint.Rdata")
load(file = "../working/dur_effects_polarization.Rdata")
load(file = "../working/dur_effects_stability.Rdata")

load(file = "../working/constraint_effects_est.Rdata")
constraint_average_treat_effects <- constraint_effects_est$average_treat_effects

load(file = "../working/polarization_effects_est.Rdata")
polarization_average_treat_effects <- polarization_effects_est$average_treat_effects

load(file = "../working/effect_estimates_stability.Rdata")
stability_average_treat_effects <- stability_effects %>% summarise(wave = "Sample One, Wave One",
                                                                   est = unique(av_est),
                                                                   high = unique(av_hi),
                                                                   low = unique(av_lo),
                                                                   name = "Trans Rights")

stability_average_treat_effects_thresh <- stability_effects_thresh %>% summarise(wave = "Sample One, Wave One",
                                                                                 est = unique(av_est),
                                                                                 high = unique(av_hi),
                                                                                 low = unique(av_lo),
                                                                                 name = "Trans Rights")




polarization_average_treat_effects_thresh$outcome <- "Polarization"
constraint_average_treat_effects_thresh$outcome <- "Constraint"
stability_average_treat_effects_thresh$outcome <- "Stability"

polarization_average_treat_effects$outcome <- "Polarization"
constraint_average_treat_effects$outcome <- "Constraint"
stability_average_treat_effects$outcome <- "Stability"

polarization_average_treat_effects_thresh <- polarization_average_treat_effects_thresh %>% filter(wave == "Sample One, Wave One")
polarization_average_treat_effects <- polarization_average_treat_effects %>% filter(wave == "Sample One, Wave One")

constraint_average_treat_effects_thresh <- constraint_average_treat_effects_thresh %>% filter(wave == "Sample One, Wave One")
constraint_average_treat_effects <- constraint_average_treat_effects %>% filter(wave == "Sample One, Wave One")

polarization_average_treat_effects_thresh$thresh <- TRUE
polarization_average_treat_effects$thresh <- FALSE
stability_average_treat_effects_thresh$thresh <- TRUE
stability_average_treat_effects$thresh <- FALSE
constraint_average_treat_effects_thresh$thresh <- TRUE
constraint_average_treat_effects$thresh <- FALSE

bind_rows(polarization_average_treat_effects_thresh,
          polarization_average_treat_effects,
          stability_average_treat_effects_thresh,
          stability_average_treat_effects,
          constraint_average_treat_effects_thresh,
          constraint_average_treat_effects) %>%
  mutate(thresh = ifelse(thresh, "> 30 seconds", "Full sample")) %>%
  ggplot(aes(x = est, xmin = low, xmax = high, y = outcome, group = thresh, col = thresh)) + 
  geom_pointrange(position = position_dodge(.5)) +
  theme_bw() +
  xlab("Average Effect of Reason-Giving \n(across all isssues)") +
  ylab("") +
  xlim(c(-.5, .5)) + 
  theme(panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        legend.background = element_rect(fill = "transparent", color = NA),
        legend.box.background = element_rect(fill = "transparent", color = NA),
        legend.key = element_blank()) + 
  scale_color_manual("Treatment group", values = c("black", "gray")) + 
  geom_vline(xintercept = 0, lty = 2)

ggsave("../out/outcomes/heterogeneous_effects_treatment_duration.pdf", width = 7, height = 4)
